Maps in R!
ST465– Making Maps in R
Today we’ll go over how to make maps in R and how to work with
spatial data sets. Some of the key R packages that we’ll work with are
tmap, mapview for map making and
sf, spData, spDataLarge for
working with spatial data.
## R packages
setwd(getwd())
library(rmarkdown)
library(tmap) ## Very commonly used -- to make static and interactive plots
library(mapview)
library(leaflet)
library(shiny)
library(tidyverse)
library(spDataLarge) ## Data sets -- use commands in notes to download
library(spData) ## Data sets
library(sf) ## simple features tmap
Let’s work with the world and coffee data
sets. For each, identify the type of object it is and its coordinate
reference system (CRS).
#data("world")
world## Simple feature collection with 177 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 11
## iso_a2 name_l~1 conti~2 regio~3 subre~4 type area_~5 pop lifeExp gdpPe~6
## * <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 FJ Fiji Oceania Oceania Melane~ Sove~ 1.93e4 8.86e5 70.0 8222.
## 2 TZ Tanzania Africa Africa Easter~ Sove~ 9.33e5 5.22e7 64.2 2402.
## 3 EH Western~ Africa Africa Northe~ Inde~ 9.63e4 NA NA NA
## 4 CA Canada North ~ Americ~ Northe~ Sove~ 1.00e7 3.55e7 82.0 43079.
## 5 US United ~ North ~ Americ~ Northe~ Coun~ 9.51e6 3.19e8 78.8 51922.
## 6 KZ Kazakhs~ Asia Asia Centra~ Sove~ 2.73e6 1.73e7 71.6 23587.
## 7 UZ Uzbekis~ Asia Asia Centra~ Sove~ 4.61e5 3.08e7 71.0 5371.
## 8 PG Papua N~ Oceania Oceania Melane~ Sove~ 4.65e5 7.76e6 65.2 3709.
## 9 ID Indones~ Asia Asia South-~ Sove~ 1.82e6 2.55e8 68.9 10003.
## 10 AR Argenti~ South ~ Americ~ South ~ Sove~ 2.78e6 4.30e7 76.3 18798.
## # ... with 167 more rows, 1 more variable: geom <MULTIPOLYGON [arc_degree]>,
## # and abbreviated variable names 1: name_long, 2: continent, 3: region_un,
## # 4: subregion, 5: area_km2, 6: gdpPercapcoffee_data## # A tibble: 47 x 3
## name_long coffee_production_2016 coffee_production_2017
## <chr> <int> <int>
## 1 "Angola" NA NA
## 2 "Bolivia" 3 4
## 3 "Brazil" 3277 2786
## 4 "Burundi" 37 38
## 5 "Cameroon" 8 6
## 6 "Central African Republic" NA NA
## 7 "Congo, Dem. Rep. of" 4 12
## 8 "Colombia" 1330 1169
## 9 "Costa Rica" 28 32
## 10 "C\u00f4te d'Ivoire" 114 130
## # ... with 37 more rowsworld_coffee = dplyr::left_join(world, coffee_data, by = "name_long")
world_coffee## Simple feature collection with 177 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## # A tibble: 177 x 13
## iso_a2 name_l~1 conti~2 regio~3 subre~4 type area_~5 pop lifeExp gdpPe~6
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 FJ Fiji Oceania Oceania Melane~ Sove~ 1.93e4 8.86e5 70.0 8222.
## 2 TZ Tanzania Africa Africa Easter~ Sove~ 9.33e5 5.22e7 64.2 2402.
## 3 EH Western~ Africa Africa Northe~ Inde~ 9.63e4 NA NA NA
## 4 CA Canada North ~ Americ~ Northe~ Sove~ 1.00e7 3.55e7 82.0 43079.
## 5 US United ~ North ~ Americ~ Northe~ Coun~ 9.51e6 3.19e8 78.8 51922.
## 6 KZ Kazakhs~ Asia Asia Centra~ Sove~ 2.73e6 1.73e7 71.6 23587.
## 7 UZ Uzbekis~ Asia Asia Centra~ Sove~ 4.61e5 3.08e7 71.0 5371.
## 8 PG Papua N~ Oceania Oceania Melane~ Sove~ 4.65e5 7.76e6 65.2 3709.
## 9 ID Indones~ Asia Asia South-~ Sove~ 1.82e6 2.55e8 68.9 10003.
## 10 AR Argenti~ South ~ Americ~ South ~ Sove~ 2.78e6 4.30e7 76.3 18798.
## # ... with 167 more rows, 3 more variables: geom <MULTIPOLYGON [arc_degree]>,
## # coffee_production_2016 <int>, coffee_production_2017 <int>, and abbreviated
## # variable names 1: name_long, 2: continent, 3: region_un, 4: subregion,
## # 5: area_km2, 6: gdpPercapWe’ve created a new data set called world_coffee, what type of object is it and what is its CRS?
Let’s plot the spatial object using the R package
tmap.
## tmap_mode "view" set to interactive viewing
tmap_mode("plot")
facets = c("coffee_production_2016", "coffee_production_2017")
tm_shape(world_coffee) +
tm_polygons(facets) +
tm_facets(nrow = 1, sync = TRUE)
Now let’s use the World data set (not that
world and World are different data sets.)
data("World")
World ## Simple feature collection with 177 features and 15 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.9 xmax: 180 ymax: 83.64513
## Geodetic CRS: WGS 84
## First 10 features:
## iso_a3 name sovereignt continent
## 1 AFG Afghanistan Afghanistan Asia
## 2 AGO Angola Angola Africa
## 3 ALB Albania Albania Europe
## 4 ARE United Arab Emirates United Arab Emirates Asia
## 5 ARG Argentina Argentina South America
## 6 ARM Armenia Armenia Asia
## 7 ATA Antarctica Antarctica Antarctica
## 8 ATF Fr. S. Antarctic Lands France Seven seas (open ocean)
## 9 AUS Australia Australia Oceania
## 10 AUT Austria Austria Europe
## area pop_est pop_est_dens economy
## 1 652860.000 [km^2] 28400000 4.350090e+01 7. Least developed region
## 2 1246700.000 [km^2] 12799293 1.026654e+01 7. Least developed region
## 3 27400.000 [km^2] 3639453 1.328268e+02 6. Developing region
## 4 71252.172 [km^2] 4798491 6.734519e+01 6. Developing region
## 5 2736690.000 [km^2] 40913584 1.495003e+01 5. Emerging region: G20
## 6 28470.000 [km^2] 2967004 1.042151e+02 6. Developing region
## 7 12259213.973 [km^2] 3802 3.101341e-04 6. Developing region
## 8 7257.455 [km^2] 140 1.929051e-02 6. Developing region
## 9 7682300.000 [km^2] 21262641 2.767744e+00 2. Developed region: nonG7
## 10 82523.000 [km^2] 8210281 9.949082e+01 2. Developed region: nonG7
## income_grp gdp_cap_est life_exp well_being footprint inequality
## 1 5. Low income 784.1549 59.668 3.8 0.79 0.42655744
## 2 3. Upper middle income 8617.6635 NA NA NA NA
## 3 4. Lower middle income 5992.6588 77.347 5.5 2.21 0.16513372
## 4 2. High income: nonOECD 38407.9078 NA NA NA NA
## 5 3. Upper middle income 14027.1261 75.927 6.5 3.14 0.16423830
## 6 4. Lower middle income 6326.2469 74.446 4.3 2.23 0.21664810
## 7 2. High income: nonOECD 200000.0000 NA NA NA NA
## 8 2. High income: nonOECD 114285.7143 NA NA NA NA
## 9 1. High income: OECD 37634.0832 82.052 7.2 9.31 0.08067825
## 10 1. High income: OECD 40132.6093 81.004 7.4 6.06 0.07129351
## HPI geometry
## 1 20.22535 MULTIPOLYGON (((61.21082 35...
## 2 NA MULTIPOLYGON (((16.32653 -5...
## 3 36.76687 MULTIPOLYGON (((20.59025 41...
## 4 NA MULTIPOLYGON (((51.57952 24...
## 5 35.19024 MULTIPOLYGON (((-65.5 -55.2...
## 6 25.66642 MULTIPOLYGON (((43.58275 41...
## 7 NA MULTIPOLYGON (((-59.57209 -...
## 8 NA MULTIPOLYGON (((68.935 -48....
## 9 21.22897 MULTIPOLYGON (((145.398 -40...
## 10 30.47822 MULTIPOLYGON (((16.97967 48...HPI stands for “Happy Planet Index” and provides a quantification of a country’s “happiness”. More information here: https://happyplanetindex.org Let’s check it out:
tmap_mode("view")
tm_shape(World) +
tm_polygons("HPI")Let’s look at the HPI along with economy type:
tmap_mode("view")
tm_shape(World) +
tm_polygons(c("HPI", "economy")) +
tm_facets(sync = TRUE, ncol = 2)We continue to make maps, this time with the metro data
set. What information is provided? What is the geometry type and what is
the CRS?
## reading in data
data(metro)
metro ## Simple feature collection with 436 features and 12 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -123.1193 ymin: -37.814 xmax: 174.7667 ymax: 60.16925
## Geodetic CRS: WGS 84
## First 10 features:
## name name_long iso_a3 pop1950 pop1960 pop1970 pop1980
## 2 Kabul Kabul AFG 170784 285352 471891 977824
## 8 Algiers El Djazair (Algiers) DZA 516450 871636 1281127 1621442
## 13 Luanda Luanda AGO 138413 219427 459225 771349
## 16 Buenos Aires Buenos Aires ARG 5097612 6597634 8104621 9422362
## 17 Cordoba Cordoba ARG 429249 605309 809794 1009521
## 25 Rosario Rosario ARG 554483 671349 816230 953491
## 32 Yerevan Yerevan ARM 341432 537759 778158 1041587
## 33 Adelaide Adelaide AUS 429277 571822 850168 971856
## 34 Brisbane Brisbane AUS 441718 602999 904777 1134833
## 37 Melbourne Melbourne AUS 1331966 1851220 2499109 2839019
## pop1990 pop2000 pop2010 pop2020 pop2030 geometry
## 2 1549320 2401109 3722320 5721697 8279607 POINT (69.17246 34.52889)
## 8 1797068 2140577 2432023 2835218 3404575 POINT (3.04197 36.7525)
## 13 1390240 2591388 4508434 6836849 10428756 POINT (13.23432 -8.83682)
## 16 10513284 12406780 14245871 15894307 16956491 POINT (-58.40037 -34.60508)
## 17 1200168 1347561 1459268 1562509 1718192 POINT (-64.18105 -31.4135)
## 25 1083819 1152387 1298073 1453814 1606993 POINT (-60.63932 -32.94682)
## 32 1174524 1111301 1065597 1023703 1057459 POINT (44.51462 40.182)
## 33 1081618 1141623 1217990 1320783 1505422 POINT (138.5986 -34.92866)
## 34 1381306 1666203 2033617 2388517 2721325 POINT (153.0281 -27.46794)
## 37 3154314 3460541 3951216 4500501 5070873 POINT (144.9633 -37.814)tmap_mode("view")
tm_basemap("Stamen.Watercolor") +
tm_shape(metro) + tm_bubbles(size = "pop2020", col = "red") +
tm_tiles("Stamen.TonerLabels")You can also play around with the style in tmap.
tmap_style("classic")
## tmap style set to "classic"
## other available styles are: "white", "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "watercolor"
tm_shape(World) +
tm_polygons("HPI", legend.title = "Happy Planet Index")We can also look through some data from New Zealand. As the others, what type of object is it and what is the CRS? How it different from others and why might that be the case?
## New Zealand
nz## Simple feature collection with 16 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1090144 ymin: 4748537 xmax: 2089533 ymax: 6191874
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
## First 10 features:
## Name Island Land_area Population Median_income Sex_ratio
## 1 Northland North 12500.561 175500 23400 0.9424532
## 2 Auckland North 4941.573 1657200 29600 0.9442858
## 3 Waikato North 23900.036 460100 27900 0.9520500
## 4 Bay of Plenty North 12071.145 299900 26200 0.9280391
## 5 Gisborne North 8385.827 48500 24400 0.9349734
## 6 Hawke's Bay North 14137.524 164000 26100 0.9238375
## 7 Taranaki North 7254.480 118000 29100 0.9569363
## 8 Manawatu-Wanganui North 22220.608 234500 25000 0.9387734
## 9 Wellington North 8048.553 513900 32700 0.9335524
## 10 West Coast South 23245.456 32400 26900 1.0139072
## geom
## 1 MULTIPOLYGON (((1745493 600...
## 2 MULTIPOLYGON (((1803822 590...
## 3 MULTIPOLYGON (((1860345 585...
## 4 MULTIPOLYGON (((2049387 583...
## 5 MULTIPOLYGON (((2024489 567...
## 6 MULTIPOLYGON (((2024489 567...
## 7 MULTIPOLYGON (((1740438 571...
## 8 MULTIPOLYGON (((1866732 566...
## 9 MULTIPOLYGON (((1881590 548...
## 10 MULTIPOLYGON (((1557042 531...legend_title = expression("Area (km"^2*")")
map_nza = tm_shape(nz) +
tm_fill(col = "Land_area", title = legend_title) + tm_borders()
map_nza + tm_style("bw")map_nza + tm_style("classic")map_nza + tm_style("cobalt")map_nza + tm_style("col_blind")mapview
Let’s use mapview now, another common R package for
constructing spatial maps. We’ll use the trails and
franconia data set. What kind of objects are they and what
is its CRS?
trails ## Simple feature collection with 543 features and 3 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: 528933.9 ymin: 5415754 xmax: 715816.7 ymax: 5597457
## Projected CRS: WGS 84 / UTM zone 32N
## First 10 features:
## FGN FKN
## 1 003756/Kunigundenweg 003756
## 2 004037/Jakobsweg (Almerswind-Coburg-Lichtenfels-Bamberg-Nuernberg) 004037
## 3 005160/Steigerwaelder Jakobsweg (Bamberg-Uffenheim) 005160
## 4 022650/Sieben-Fluesse-Wanderweg 022650
## 5 005316/Rennsteigverein 1896 e.V. / Bamberger Rennsteig 005316
## 6 012029/Steigerwald-Panoramaweg 012029
## 7 013633/Frankenwaldverein / Markgrafenweg 013633
## 8 023964/Rot-Main-Auen-Weg 023964
## 9 007443/Jean-Paul-Weg 007443
## 10 007443/Jean-Paul-Weg 007443
## district geometry
## 1 Oberfranken MULTILINESTRING ((631670.9 ...
## 2 Oberfranken MULTILINESTRING ((636557.1 ...
## 3 Oberfranken MULTILINESTRING ((633443.8 ...
## 4 Oberfranken MULTILINESTRING ((631625.6 ...
## 5 Oberfranken MULTILINESTRING ((636894.6 ...
## 6 Oberfranken MULTILINESTRING ((634029.7 ...
## 7 Oberfranken MULTILINESTRING ((684491.3 ...
## 8 Oberfranken MULTILINESTRING ((686023 55...
## 9 Oberfranken MULTILINESTRING ((684861.6 ...
## 10 Oberfranken MULTILINESTRING ((684905.6 ...franconia## Simple feature collection with 37 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 8.975926 ymin: 48.8625 xmax: 12.27535 ymax: 50.56422
## Geodetic CRS: WGS 84
## First 10 features:
## NUTS_ID SHAPE_AREA SHAPE_LEN CNTR_CODE NAME_ASCI
## 1 DE241 0.006736012 0.3926225 DE Bamberg, Kreisfreie Stadt
## 2 DE242 0.008424469 0.6247263 DE Bayreuth, Kreisfreie Stadt
## 3 DE243 0.005982341 0.5185471 DE Coburg, Kreisfreie Stadt
## 4 DE244 0.007329480 0.4569815 DE Hof, Kreisfreie Stadt
## 5 DE245 0.146698316 3.4819699 DE Bamberg, Landkreis
## 6 DE246 0.159489736 3.6242023 DE Bayreuth, Landkreis
## 7 DE247 0.074698748 2.6954234 DE Coburg, Landkreis
## 8 DE248 0.079746707 1.7712298 DE Forchheim
## 9 DE249 0.112934151 2.7544708 DE Hof, Landkreis
## 10 DE24A 0.081960299 1.9393830 DE Kronach
## geometry district
## 1 MULTIPOLYGON (((10.92581 49... Oberfranken
## 2 MULTIPOLYGON (((11.58157 49... Oberfranken
## 3 MULTIPOLYGON (((10.95355 50... Oberfranken
## 4 MULTIPOLYGON (((11.93067 50... Oberfranken
## 5 MULTIPOLYGON (((10.87615 50... Oberfranken
## 6 MULTIPOLYGON (((11.70657 50... Oberfranken
## 7 MULTIPOLYGON (((10.88654 50... Oberfranken
## 8 MULTIPOLYGON (((11.26376 49... Oberfranken
## 9 MULTIPOLYGON (((11.91988 50... Oberfranken
## 10 MULTIPOLYGON (((11.36979 50... OberfrankenWhat does it look like if we only plot the trails data?
trails%>%
mapview()What about the franconia data set?
franconia%>%
mapview()We’d like to plot the trails onto the map, but they’re not the same
CRS. In lecture, we saw that we can use a suite of functions starting
with st_ to modify a spatial object. We’ll use a few here
to change the crs and plot the data sets together:
## first, using st_transform() function we changed the CRS
trails %>%
st_transform(st_crs(franconia)) %>%
st_intersection(franconia[franconia$district == "Oberfranken", ]) %>%
st_collection_extract("LINE") %>%
mapview(color = "red", lwd = 3, layer.name = "trails") +
mapview(franconia, zcol = "district", burst = TRUE) Lastly, suppose you also wanted to know about any breweries along the
trails. Let’s check out the breweries data set. What kind
of object is it and what is its CRS?
breweriesWhat does a map of the breweries look like?
breweries%>%
mapview()Let’s combine them all:
trails %>%
st_transform(st_crs(franconia)) %>%
st_intersection(franconia[franconia$district == "Oberfranken", ]) %>%
st_collection_extract("LINE") %>%
mapview(color = "red", lwd = 3, layer.name = "trails") +
mapview(franconia, zcol = "district", burst = TRUE) +
breweriestmap
What if we tried to do this in tmap?
tm_shape(franconia) + tm_fill("district") +
tm_shape(trails) + tm_lines() +
tm_shape(breweries) + tm_dots()